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LONG-TERM  GOALS 

Estimation  of  seabed  geoacoustic  parameters  in  shallow  water  by  acoustic  remote  sensing  remains  a 
challenging  task  due  to  constraints  on  hardware,  data  collection  and  analysis,  and  cost  of  maritime 
surveys.  This  work  focuses  on  the  application  of  two  techniques  that  might  offer  a  solution  to  those 
constraints:  the  use  of  ambient  noise  to  probe  the  seabed,  and  Bayesian  inversion  of  these  data  to 
estimate  geoacoustic  parameters  of  interest  together  with  their  uncertainties.  The  long-term  goal  of  this 
work  is  to  establish  general  methods  for  processing  and  inverting  ambient  noise  data  and  assessing  the 
quality  of  the  results  by  quantifying  their  uncertainties. 
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OBJECTIVES 


This  work  has  three  main  objectives:  First,  quantifying  the  ability  to  resolve  seabed  geoacoustic 
parameters  using  ambient  noise  measurements.  Second,  comparing  those  estimates  to  the  ones 
obtained  from  active  source  inversion  methods.  Third,  increasing  the  understanding  of  the 
experimental  conditions  and  equipment  required  for  the  collection  of  ambient  noise  data  suitable  for 
geoacoustic  inversion. 

APPROACH 

Traditional  investigation  of  seabed  sediment  properties  has  relied  heavily  on  direct  measurements,  such 
as  core  sampling  and  geo-probes,  or  indirect  measurements  with  active  systems.  Direct  methods  have 
the  evident  problem  of  lack  of  spatial  resolution  due  to  time/cost  constraints,  while  active  methods  can 
be  limited  due  to  deployment  procedures  and  environmental  concerns,  often  requiring  the  use  of  a 
vessel  to  tow  the  active  device  over  a  geographic  area  of  interest.  As  alternative  to  active  systems,  it  is 
known  that  the  wind-driven  ambient  noise  field  recorded  at  a  vertical  array  carries  infonnation  of  the 
seabed  layering  structure,1  which  is  exploited  in  this  research.  The  approach  for  this  work  consists  of: 

1)  Defining  an  inversion  method:  the  Bayesian  framework  has  been  selected  in  this  study  to  carry  out 
geoacoustic  inversion  from  ambient  noise  data.  The  work  by  Dettmer  and  Dosso3'4  (University  of 
Victoria)  and  Holland  (Pennsylvania  State  University)  on  Bayesian  controlled-source  reflection 
coefficient  inversion  is  directly  applicable  to  the  proposed  inversion  of  similar  data  as  extracted  from 
the  ambient  noise  field.5 

2)  Implementing  algorithms  for  numerical  estimation  of  the  posterior  probability  density  (PPD)  over 
the  geoacoustic  parameters  of  interest:  Since  analytical  solutions  for  the  PPD  are  generally  not 
available  for  non-linear  problems,  Markov  chain  Monte  Carlo  (MCMC)  methods  are  used  to  sample 
from  this  distribution.4  In  this  work,  Metropolis-Hastings  sampling  (MHS)  is  applied  to  detennine 
marginal  probability  densities.  Perturbations  are  applied  in  a  principal-component  parameter  space, 
which  is  a  rotated  representation  of  the  physical  parameter  space  in  which  the  axes  align  with  the 
dominant  correlation  directions.  This  rotation  provides  a  more  efficient  exploration  of  the  parameter 
space,  and  is  particularly  effective  when  strong  correlations  between  parameters  are  present. 

3)  Identifying  a  forward  model:  the  input  to  the  Bayesian  inversion  is  the  seabed  power  reflection 
coefficient  R  or  the  bottom  loss  (BL  =  -lOlogk),  which  can  be  computed  from  the  ratio  of  upward  to 
downward  energy  fluxes  obtained  by  beamforming  ambient  noise  measured  at  a  vertical  linear  array 
(VLA).2  6  The  forward  model  consists  of  computing  a  representation  of  the  ambient  noise  data 
covariance  matrix,  from  which  replicas  of  the  BL  can  be  calculated  for  different  combinations  of  the 
geoacoustic  parameters.  This  replica  BL  is  adjusted  to  include  the  smearing  effect  introduced  by  the 
VLA's  finite  aperture.  Software  routines  for  the  forward  model  have  been  validated  by  comparison 
with  OASN,  the  ambient  noise  module  from  the  wavenumber-integration  model  OASES  (OASN  itself 
is  too  computationally  expensive  to  use  in  the  inversion  algorithm). 

4)  Detennining  the  impact  of  array  design  (e.g.,  aperture,  sensor  density)  and  experimental  conditions 
(e.g.,  wind  speed)  in  geoacoustic  inversion  by  analysis  of  synthetic  data  obtained  from  the  forward 
model.5 
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5)  Applying  the  inversion  framework  to  experimental  data  from  previous  publications.  ’  ’ 
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WORK  COMPLETED 


2 

1)  The  ray-tracing  representation  of  the  ambient  noise  field  developed  by  Harrison  has  been 
adopted  as  the  forward  model  to  compute  the  angle-  and  frequency-dependent  seabed  BL.  This 
approach  considers  wind-driven  surface  dipoles  as  the  driving  mechanism  for  the  ambient  noise. 
The  strength  of  this  field  relative  to  other  unwanted  noise  mechanisms  defines  a  signal-to-noise 
ratio  (SNR),5  which  is  included  in  this  work  as  an  unknown  frequency-dependent  parameter. 

2)  Trans-dimensional  (trans-D)  Bayesian  inversion  with  parallel  tempering10  was  used  for 
geoacoustic  parameter  estimation  from  BL  data  derived  from  simulated  ambient-noise. 1 1 

3)  The  trans-D  inversion  was  applied  to  data  from  the  MAPEX  2000  experiment1 1  and  the  results 
were  compared  to  previous  work  that  utilized  the  Bayesian  information  criterion  (BIC)  for  fixed¬ 
dimensional  inversion.  Preliminary  inversions  using  data  from  a  drifting  array  from  the  Boundary 
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2003  experiment  were  also  carried  out. 
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4)  A  sequential  trans-D  Monte  Carlo  algorithm  ”  was  applied  to  simulated  data  corresponding  to  a 
drifting  vertical  array.  This  approach  provides  estimates  of  geoacoustic  parameters,  true-depth 
layering  structure  of  the  seabed,  and  parameter  uncertainties.  The  inversion  was  applied  to 
incoherent  estimates  of  seabed  power  reflection  coefficient  data,  computed  as  the  array  drifts 
along  a  range-dependent  track. 

5)  Using  simulated  data,  the  impact  of  array  aperture  in  geoacoustic  resolution  was  studied.  To 
resolve  complicated  seabed  structure,  large  apertures  are  required.  Since  minimizing  array 
aperture  offers  advantages  in  array  design  and  deployment,  signal  processing  techniques  to  extend 
short  apertures  (i.e.,  synthetic  array  apertures)  were  evaluated. 

RESULTS 

Trans-D  approach  to  model  selection:  Initial  results  in  this  effort3  illustrated  the  application  of 
Bayesian  inversion  to  BL  data  using  the  BIC  for  model  selection.  The  impact  of  wind  strength  in  the 
estimation  of  seabed  geoacoustic  parameters  was  quantified,  and  results  from  geoacoustic  inversion 
were  compared  to  direct  (core)  measurements.  As  a  more  general  approach,  the  trans-D  Bayesian 
inversion  algorithm1 1-13  was  used  to  provide  automated  model  selection  over  an  unknown  number  of 
seabed  layers  and  to  quantify  the  uncertainty  due  to  model  selection.  With  the  trans-D  method,  models 
from  a  set  of  K  candidates  are  included  in  the  estimation  of  the  PPD,  defined  as 

H (1) 

v  P(d) 

where  Z.(d  |  m k,Ik)  is  the  likelihood  function,  while P(/t)is  the  prior  distribution  for  parametrization  I  k , 
assumed  here  as  a  discrete  unifonn  distribution.  p(mk  \  Ik )  is  the  prior  distribution  for  the  geoacoustic 
parameters  m/cfor  a  layered  seabed  with  k  interfaces.  The  vector  is  defined  as 

m k  =  [c,  px  a,  z, ... ck+]  pk+]  ak+]  SNR , . . . SNRF f ,  (2) 

where  c;,  pn  an  and  z,  are  the  sound  speed,  density,  attenuation  and  interface  depth  of  the  ith  layer, 
respectively.  The  SNRs5  account  for  the  unknown  strength  of  the  wind-driven  ambient -noise  data  (i.e. 
the  useful  signal)  versus  other  unwanted  sources  of  noise  at  F  frequencies.  The  PPD  is  sampled  by  a 
reversible-jump  Markov  chain  Monte  Carlo  (rjMCMC)  algorithm,414  which  uses  an  extended 
Metropolis-Hasting  (MH)  criterion  that  allows  trans-D  jumps  between  parameterizations  4,  sampling 
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probabilistically  over  models  with  different  parametrizations  (number  of  layers)  and  quantifying  the 
uncertainty  due  to  the  lack  of  knowledge  of  the  model  parameterization. 

Inversion  results  using  a  simulated  environment  consisting  of  depth-dependent  variations  in  sound 
speed,  density,  and  attenuation11  were  in  good  agreement  with  the  true  sediment  profdes,  as  well  as 
with  similar  inversions  carried  out  with  active  source  simulated  data.  The  trans-D  inversion  was  also 
applied  to  data  obtained  from  a  moored  VLA  during  the  MAPEX  2000  experiment,11  providing 
improved  estimates  of  parameter  uncertainties  compared  to  previous  work  under  this  effort5  using  the 
BIC  approach  to  model  selection. 


Inversion  of  drifting-array  simulated  data:  For  a  drifting  array,  the  PPD  evolves  with  time  as  the 
array  moves  over  sediments  in  which  the  number  of  layers,  the  depth  of  interfaces,  or  the  geoacoustic 
parameters  change  as  a  function  of  range.  Sequential  datasets  can  then  be  obtained  by  discretizing 
continuous-time  recordings  of  ambient  noise.  For  this  application,  a  particle  filter  ~  is  used  to  update 
the  estimated  geoacoustic  parameters  from  one  array  position  to  the  next  as  new  data  become 
available.  To  generate  simulated  sequential  data,  the  environment  shown  in  Fig.  l-(a)  (similar  panels  for 
density  and  sediment  attenuation  were  included  in  the  simulation)  was  input  to  OASES7  for 
computation  of  the  range-dependent  ambient-noise  field  at  a  32-element  VLA  with  0.18  m  inter- 
element  spacing.  This  simulated  enviromnent  (used  in  previous  work  in  the  context  of  active -source 
surveys  )  includes  realistic  features  such  as  range-dependent  smooth  transitions  in  geoacoustic 
parameters,  thin  sediment  layers,  and  abrupt  variations  introduced  by  a  geologic  fault  and  an  erosional 
channel.  Conventional  beamfonning  was  used  to  obtain  the  power  reflection  coefficient  at  8 
frequencies  from  1000  Hz  to  3000  Hz  and  grazing  angles  from  20°  -  90°.  These  data  were  provided  to 
the  sequential  Bayesian  trans-D  Monte  Carlo  algorithm  for  estimation  of  the  PPD.  Figure  l-(a) 
(bottom)  shows  the  estimated  mean  model  of  the  sediment  sound  speed,  corresponding  to  simulated 
data  from  a  224-element  array  (40  m  aperture). 


(a) 


(b)  Power  reflection  coefficient  (150  m  range) 


True  seabed  sound  speed  profile 


Mean  model  via  inversion 


20  40  60  80  100  120  140  160 

Range (m) 


Sound 

speed 

(m/s) 

1 1700 

1600 

y  isoo 


N 

I 


N 

X 


u 

3 

= 

O- 

o 


[i. 


2.5 
2.0 

1.5 
1 

From  mean  model 


From  true  model 


Magnitude 

|f 

0.5 

i 


0  20  40  60  80 

Grazing  angle  (degrees) 


Figure  1  Geoacoustic  inversion  of  ambient-noise  sequential  data:  (a)  True  seabed  sound  speed 
profile  input  to  OASES  to  generate  simulated  data  vs  mean  model  obtained  via  inversion,  (b)  Power 
reflection  coefficient  estimated  by  beamforming  using  the  true  parameters  (top)  from  (a)  and  the 
most  probable  model  obtained  via  inversion  (bottom),  at  range  150  m. 
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The  estimated  geoacoustic  parameters  (sound  speed,  number/depth  of  sediment  interfaces)  closely 
resemble  the  true  sediment  profile  even  at  ranges  containing  thin  sediment  layers  which  are 
challenging  for  an  inversion  algorithm.  Similar  agreement  was  obtained  for  the  density  and  the 
attenuation. 


As  an  example  of  the  resolution  limitation  of  the  data,  Fig.  1  -(a)  exhibits  a  mismatch  in  the  inversion 
results  at  ranges  148-154  m  corresponding  to  the  erosional  channel.  To  explain  this  mismatch, 

Fig.  1  -(b)  shows  the  power  reflection  coefficient  estimated  by  beamforming  using  the  true  geoacoustic 
parameters  (top)  and  the  most  likely  model  obtained  via  inversion  (bottom).  The  structure  of  the 
interference  pattern  in  both  cases  is  very  similar  for  the  angular  range  used  for  inversion  (20°  -  90°), 
causing  the  inversion  algorithm  to  have  slow  convergence  to  the  true  solution  while  searching  over  the 
parameter  space.  Given  more  iterations,  the  estimated  profile  should  eventually  converge  to  the 
theoretical  environment.  Figure  1  shows  the  best  results  in  a  series  of  simulations  in  which  the  length 
of  the  array  was  varied  from  40  m  to  6  m  to  observe  the  impact  of  array  aperture  into  the  estimated 
sediment  profiles.  As  the  aperture  decreased,  features  in  the  power  reflection  coefficient  are  lost  due  to 
beamforming  smearing,  resulting  in  low  geoacoustic  resolution.  For  field  experiments,  increasing 
geoacoustic  resolution  by  augmenting  the  array  aperture  would  make  array  deployment  more  difficult 
and  decrease  array  stability  while  drifting.  Therefore,  other  options  such  as  synthetically  extending  the 
array  were  analyzed. 


Extension  of  the  array  aperture  by  signal  processing:  Given  an  A-element  array  with  inter-element 
spacing  Az ,  the  seabed  power  reflection  coefficient  at  grazing  angle  0  can  be  approximated 
as  R(6)  =  B(-0)/B(6),  where  B(0)  is  the  output  of  the  conventional  beamformer: 


B(G)  =  wH(0) 
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w($)  =  wH(6)  Cw(6*).  (3) 


In  (3),  w(0)  is  the  array  steering  vector,  while  Cba  is  the  spatial  coherence  between  two  sensors 
separated  a  distance ( a  -  b)Az ,  with  /  <  b  <  a  <  N .  Under  mild  requirements15,  the  matrix  C  for 
ambient  noise  at  a  vertical  array  can  be  modelled  as  Toeplitz  and  only  the  first  row  must  be  known  to 
uniquely  determine  the  matrix.  Increasing  the  array  aperture  by  signal  processing  techniques  (as 
opposed  to  by  physically  adding  sensors  to  the  array)  requires  extension  of  the  coherence  function, 
which  can  be  achieved  by  zero  padding15  or  by  extrapolation.  For  example,  to  obtain  a  synthetic  6- 
element  array  from  a  3-element  array  it  is  required  to  approximate  the  true  spatial  coherence  as: 


[c„  cn  cu  cu  c„ 


cn  C„  0  0  0\ 

Cn  Cn  C„  C,s 


zero  -  padding, 


extrapolation. 


(4) 


The  zero-padding  method  has  been  shown  to  enhance  features  of  the  estimated  power  reflection 
coefficient.15  However,  the  power  reflection  coefficient  obtained  by  a  zero-padded  covariance  matrix 
differs  from  the  corresponding  quantity  that  would  be  obtained  by  using  an  actual  array  of  similar 
aperture.  In  inversion,  this  represents  a  potential  problem,  since  the  forward  model  used  while 
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navigating  the  parameter  space  must  accurately  represent  the  data  to  obtain  unbiased  estimation  of  the 
geoacoustic  parameters  of  interest.  To  avoid  this  issue,  for  inversion  problems  extrapolation  of  the 
coherence  function  should  be  considered. 


The  spatial  coherence  due  to  surface-generated  ambient  noise  can  be  written  as15 

Ml 

C(z)  *  2n  |  [l  -  /?(*)]  “ \eikz  +  R{k)e~ikz  }dk,  (5) 

o 

where  X  is  the  wavelength,  R(k)  is  the  seabed  power  reflection  coefficient,  and  k  =  2xsm9 1  A.  is  the 
vertical  spatial  wavenumber.  It  has  also  been  shown15  that  the  Fourier  transfonn  of  the  spatial 
coherence  is 


F{C(z)}« 


|[l -*(*)] -'*(*) 

IMWF 


-\/A<k<0, 
0<k<  -1/A. 


(6) 


Therefore,  C(z)  is  band-limited  and  methods  such  as  the  iterative  Papoulis-Gerchberg  algorithm16  or  its 
non-iterative  versions  can  be  used  for  extrapolation.  Figure  2(a)  shows  an  example  of  the  estimation 
of  the  bottom  loss  for  a  two-layer  environment  consisting  of  a  0.75  m  sediment  layer  (sound 
speed=1550  m/s,  density=1.5  g/cm attenuation=0.2  dB/7)  on  top  of  a  halfspace  (1600  m/s,  2.0  g/cm3, 
0.15  dBA). 


(a)  (b) 


Figure  2  Comparison  of  array-extension  techniques  for  a  1.5-m  aperture  array  used  to  approximate 
a  6-m  array  in  a  2-layer  environment:  (a)  Bottom  loss  obtained  from  the  analytical  solution  and 
from  beamforming  using  true  apertures  of  1.5  m  and  6  m,  as  well  as  6  m  synthetic  apertures  using 
the  zero-padding15 and  extrapolation16 approaches,  (b)  Spatial  coherence  (4)  for  a  true  6-m  aperture, 
a  zero-padded  6-m  aperture,  and  an  extrapolated  6-m  aperture. 
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Figure  2-(a)  shows  an  example  in  which  the  actual  array  is  only  1.5  m  aperture.  The  true  bottom  loss  is 
shown  as  a  reference,  with  sharp  peaks  at  18°,  30°,  and  49°.  In  this  case,  the  zero-padding  technique 
fails  to  reveal  the  BL  peaks  at  30°  and  49°  when  used  to  synthesize  a  6-m  aperture,  while  the 
extrapolation  technique  shows  those  peaks.  For  inversion  purposes,  the  most  important  characteristic 
of  the  extrapolation  technique  is  that  the  results  from  the  extrapolated  6-m  aperture  are  in  reasonable 
agreement  with  the  results  from  a  true  6-m  aperture.  Figure  2(b)  compares  the  true  coherence  function 
over  6  m,  to  the  extrapolated  coherence  function  and  the  zero-padded  coherence.  The  extrapolated 
coherence  resembles  the  true  coherence,  following  the  same  oscillatory  pattern.  The  difference 
between  the  extrapolated  and  the  true  coherence  is  due  to  the  slow  convergence  of  the  Papoulis- 
Gerchberg  algorithm,16  which  in  Figure  2(b)  was  truncated  at  1000  iterations  (~2  minutes 
computational  time). 

Figure  3 -(a)  shows  a  similar  example  with  a  true  6-m  aperture  array  (dashed  curve),  in  which  the  BL 
peaks  at  30°  and  49°  can  be  observed.  The  zero-padding  technique  was  used  to  double  the  length  of 
the  array,  and  the  result  shows  sharper  peaks  approaching  the  true  bottom  loss.  However,  the  BL  is 
underestimated  at  grazing  angles  around  90°,  and  in  general  the  results  from  this  extrapolation  do  not 
agree  with  the  results  from  a  true  12-m  aperture.  The  extrapolation  technique  using  the  Papoulis- 
Gerchberg  algorithm  is  also  demonstrated,  showing  also  an  improvement  in  the  height  of  the  BL  while 
being  in  good  agreement  with  the  true  12-m  aperture  over  most  grazing  angles. 


(a) 


(b) 


Figure  3  Comparison  of  array-extension  techniques,  similar  to  Fig.2,  with  a  6-m  aperture  array 
used  to  approximate  a  12-m  aperture  array.  As  in  Fig.2,  the  results  from  the  extrapolated  12-m 
aperture  are  in  good  agreement  with  the  results  that  would  be  obtained  with  a  true  12-m  aperture. 
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IMPACT/APPLICATIONS 


In  shallow  water  regions  the  performance  of  Navy  sonar  systems  is  strongly  influenced  by  acoustic 
interaction  with  the  seabed  and,  therefore,  knowledge  of  geoacoustic  parameters  and  their 
corresponding  uncertainties  are  required  to  predict  and  optimize  sonar  perfonnance.  Bayesian 
inversion  methods  offer  an  elegant  and  powerful  framework  not  only  for  parameter  extraction  but  also 
for  uncertainty  estimation,  thereby  quantifying  the  geoacoustic  infonnation  content  of  the  data.  The 
proposed  inversion  methodology  has  been  highly  effective  when  applied  to  active  surveys,  and  current 
results5’11'13  using  experimental  and  simulated  ambient  noise  data  show  great  potential  to  overcome 
limitations  of  current  methods  of  geoacoustic  inversion.  The  smearing  introduced  by  beamforming  is 
one  of  the  main  factors  limiting  the  resolution  of  ambient  noise  remote  sensing  methods.  Since 
increasing  the  array  length  is  an  impractical  solution  to  these  problems,  zero-padding  and  extrapolation 
techniques  may  be  helpful. 

RELATED  PROJECTS 

1)  Automated  geoacoustic  inversion  and  uncertainty:  Meso-scale  seabed  variability  in  shallow  water 
environments  (Award  Number:  N00014-09-1-0394).  This  project  develops  a  Bayesian 
methodology  for  advanced  and  automated  geoacoustic  inversion.  A  range  of  active  source  data 
are  inverted  to  quantify  geoacoustic  uncertainty.  This  project  applies  and  further  develops  these 
methods  for  ambient -noise  data. 

2)  Ocean  Ambient  Noise  Studies  for  Shallow  and  Deep  Water  Environments,  2012-2014  (Award 
Number:  N0001 4-1 2-1 -001 7). 
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